%111111111111  11111111111111111111111111111111  11111111111111111111
%111111111111  11111111111111111111111111111111  11111111111111111111
%111111111111  11111111111111111111111111111111  11111111111111111111
%111111111111  11111111111111111111111111111111  11111111111111111111
%111111111        11111111111111111111111111        11111111111111111
%111111111   b->  11111111111111111111111111        11111111111111111
%              a---->
%111111111        11111111111111111111111111        11111111111111111
%111111111        11111111111111111111111111        11111111111111111
%111111111111  11111111111111111111111111111111  1111111111111
%111111111111  11111111111111111111111111111111  1111111111111
%111111111111  11111111111111111111111111111111  1111111111111
%111111111111  1<-----------c---------------->1  1111111111111
%111111111111  11111111111111111111111111111111  1111111111111
%111111111111  11111111111111111111111111111111  1111111111111
%111111111111  11111111111111111111111111111111  1111111111111   depth=d

clear;



d=10;
b=d/5*8;
a=d/5*24;
c=d/5*50;
left=round(c*1);    %=========parameter==========
right=round(c*1);
up=round(c*0.8);
down=round(c*1);
reservior=120;
lx=4;
ly=5;
%======Random==channels=====
ch1=d/5*3;
ch2=d/5*5;
ch3=d/5*6.4;
ch4=d/5*10.4;

%ch_int=[round(ch1/2),round(ch2/2),round(ch3/2),round(ch4/2)];
%ch_int=[floor(ch1/2),floor(ch2/2),floor(ch3/2),floor(ch4/2)];
ch_int=[ceil(ch1/2),ceil(ch2/2),ceil(ch3/2),ceil(ch4/2)];
%===========================
% a2=round(a/2);
% b2=round(b/2);
% c2=round(c/2);


a2=floor(a/2);
b2=floor(b/2);
c2=floor(c/2);
%===============================



%===========RANDOM PORE SIZE=============
pore_size=[18,24,36];%radius of pore
[rand_x,nrand]=size(pore_size);

%========================================





xdir=[left];
ydir=[down];

for i=1:lx-1
  xdir=[xdir,c];
end
for i=1:ly-1
  ydir=[ydir,c];
end
xdir=[xdir,right];
ydir=[ydir,up];

%ydir=[down,c,c,c,up];
%xdir=[left,c,c,c,right];
nx=left+c*(lx-1)+right;
ny=up+down+reservior+c*(ly-1);

pa=ones(nx,ny);


for i=1:nx
    for j=1:reservior
        pa(i,j)=0;
    end
end

%===========make pores==================
for indx=1:lx
    for indy=1:ly
        corx=(indx-1)*c+left;
        cory=reservior+(indy-1)*c+down;
        %=======random pore size======
        a2=pore_size(1,fix(rand(1)*nrand+1));
        %=============================
        for (i=corx-a2:corx+a2)
            for (j=cory-a2:cory+a2)
                pa(i,j)=0;
            end
        end
        
    
    end
end
%=========================================

%===============equal===channels==========
% for indx=1:lx
%     indy=1
%         corx=(indx-1)*c+left;
%         cory=reservior+(indy-1)*c+down;
%         for (i=corx-b2:corx+b2)
%             for j=1:ny
%                 pa(i,j)=0;
%             end
%         end 
% end
%     
% for indy=1:ly
%     indx=1
%         corx=(indx-1)*c+left;
%         cory=reservior+(indy-1)*c+down;
%         for (j=cory-b2:cory+b2)
%             for i=1:nx
%                 pa(i,j)=0;
%             end
%         end 
% end

%============================================



%===============Random===channels==========
texts=[];
cory2=reservior;
cory1=reservior;
for indy=1:ly+1
    cory2=cory2+ydir(indy);
    
    for indx=1:lx
    
        corx=(indx-1)*c+left;
    
        chw=ceil(rand*4);
        for i=corx-ch_int(chw):corx+ch_int(chw)
            for j=cory1:cory2
                pa(i,j)=0;
            end
           
        end
        texts=[texts;[corx-ch_int(chw)-10,cory1+abs(cory1-cory2)/2,2*ch_int(chw)]];
    end
    cory1=cory1+ydir(indy);
end


corx2=1;
corx1=1;
for indx=1:lx+1
    corx2=corx2+xdir(indx);
    
    for indy=1:ly
    
        cory=reservior+(indy-1)*c+down;
    
        chw=ceil(rand*4);
        for i=corx1:corx2
            for j=cory-ch_int(chw):cory+ch_int(chw)
                pa(i,j)=0;
            end
           
        end
        texts=[texts;[corx1+abs(corx1-corx2)/2,cory+ch_int(chw)+10,2*ch_int(chw)]];
    end
    corx1=corx1+xdir(indx);
end
%============================================


%===============side channels================
% side_w=9;
% for i=1:side_w
%     for j=reservior+down:ny
%         pa(i,j)=0;
%     end
% end
% 
% for i=nx-side_w+1:nx
%     for j=reservior+down:ny
%         pa(i,j)=0;
%     end
% end
%============================================

imshow(pa);

figure
imshow(pa);
[tx,ty]=size(texts);
for i=1:tx
    text(texts(i,2),texts(i,1),int2str(texts(i,3)),'FontSize',13);
end
%============================================
% fid = fopen('Emily_JunctionS.dat','wt');
% 
% for k=1:d+2
%     for j=1:ny
%         for i=1:nx
%         if ((k==1) || (k==d+2))
%             if (j<=reservior)
%                 fprintf(fid,'%1d\n',0);
%             else
%                 fprintf(fid,'%1d\n',1);
%             end
%         else
%             fprintf(fid,'%1d\n',pa(i,j));
%         end
%         end
%     end
% end
%=============================================

nx
ny
d+2
